Drought prediction in continental Italy: an SpVAR approach
We set the seed for reproducibility and import the necessary libraries and functions.
set.seed(123)
suppressPackageStartupMessages({
library(dplyr)
library(tidyr)
library(ggplot2)
library(systemfit)
library(plm)
library(ncdf4)
library(abind)
library(ggmap)
library(Rmpfr)
library(gridExtra)
library(geometry)
library(caret)
library(cowplot)
library(forecast)
})
setwd(
"/Users/enrico/Library/Mobile Documents/com~apple~CloudDocs/Documents/Carriera/Scuola:Università/Bocconi Master/Thesis"
)
source("functions.r")We read the SPEI data and prepare it.
# Open the netCDF file with monthly SPEI data
SPEI <- nc_open("Data/spei12.nc")
# Get the data from the first variable
SPEI_data <- ncvar_get(SPEI, attributes(SPEI$var)$names[1])
# Get the longitude and latitude values
SPEI_lon <- ncvar_get(SPEI, attributes(SPEI$dim)$names[1])
SPEI_lat <- ncvar_get(SPEI, attributes(SPEI$dim)$names[2])
# Close the file
nc_close(SPEI)
# Subset the data for the relevant time steps
SPEI_last_year<-SPEI_data[,,1367:1368]
SPEI_relevant_years <- SPEI_data[, , 967:1368]We read the SMoist data and prepare it.
# Open the netCDF file with Soil Moisture data
SMoist <- nc_open("Data/SoilMoisture.nc")
# Get the data from the first variable
SMoist_data <- ncvar_get(SMoist, attributes(SMoist$var)$names[1])
# Order the data with the same ordering as the SPEI
SMoist_ordered <- SMoist_data[361:720, ,]
SMoist_ordered <-
abind(SMoist_ordered, SMoist_data[1:360, ,], along = 1)
# Get the longitude and latitude values
SMoist_lon <- ncvar_get(SMoist, attributes(SMoist$dim)$names[1])
SMoist_lat <- ncvar_get(SMoist, attributes(SMoist$dim)$names[3])
# Print a message to check if the longitudes and latitudes match with the ones from SPEI
print(paste(
"Check on longitudes",
toString(all(SMoist_lon - 180 == SPEI_lon)),
"check on latitudes",
toString(all(SMoist_lat == SPEI_lat))
))## [1] "Check on longitudes TRUE check on latitudes TRUE"# Close the file
nc_close(SMoist)
# Subset the data for the relevant time steps
SMoist_last_year <- SMoist_ordered[, , 659:660]
SMoist_relevant_years <- SMoist_ordered[, , 259:660]We read the NDVI data and prepare it.
# Open the netCDF file with NDVI data
NDVI_read <- nc_open("Data/NDVI.nc")
# Get the data from the first variable
NDVI <- ncvar_get(NDVI_read, attributes(NDVI_read$var)$names[1])
# Get the longitude and latitude values
NDVI_lon <- ncvar_get(NDVI_read, attributes(NDVI_read$dim)$names[1])
NDVI_lat <- ncvar_get(NDVI_read, attributes(NDVI_read$dim)$names[2])
# Print a message to check if the longitudes and latitudes match with the ones from SPEI
print(paste(
"Check on longitudes",
all(NDVI_lon == SPEI_lon),
"check on latitudes",
all(NDVI_lat == SPEI_lat)
))## [1] "Check on longitudes TRUE check on latitudes TRUE"# Close the file
nc_close(NDVI_read)
# Subset the data for the relevant time steps
NDVI_last_year <- NDVI[, , 401:402]We read the PET data and prepare it.
# Open the first netCDF file with monthly PET data and store it
PET <- nc_open("Data/cru_ts4.03.1981.1990.pet.dat.nc")
PET_data <- ncvar_get(PET, attributes(PET$var)$names[1])
# Subset the data for the relevant time steps
PET_data <- PET_data[, , 7:120]
# Close the file
nc_close(PET)
# Open the second netCDF file with monthly PET data and store it
PET <- nc_open("Data/cru_ts4.03.1991.2000.pet.dat.nc")
PET_data_2 <- ncvar_get(PET, attributes(PET$var)$names[1])
# Close the file
nc_close(PET)
# Open the third netCDF file with monthly PET data and store it
PET <- nc_open("Data/cru_ts4.03.2001.2010.pet.dat.nc")
PET_data_3 <- ncvar_get(PET, attributes(PET$var)$names[1])
# Close the file
nc_close(PET)
# Open the fourth netCDF file with monthly PET data and store it
PET <- nc_open("Data/cru_ts4.03.2011.2018.pet.dat.nc")
PET_data_4 <- ncvar_get(PET, attributes(PET$var)$names[1])
# Subset the data for the relevant time steps
PET_data_4 <- PET_data_4[, , 1:48]
# Get the longitude and latitude values
PET_lon <- ncvar_get(PET, attributes(PET$dim)$names[1])
PET_lat <- ncvar_get(PET, attributes(PET$dim)$names[2])
# Print a message to check if the longitudes and latitudes match with the ones from SPEI
print(paste(
"Check on longitudes",
all(PET_lon == SPEI_lon),
"check on latitudes",
all(PET_lat == SPEI_lat)
))## [1] "Check on longitudes TRUE check on latitudes TRUE"# Close the file
nc_close(PET)
# Merge the 4 PET datasets together
PET <- array(0, dim = c(720, 360, 402))
PET[, , ] <- c(PET_data, PET_data_2, PET_data_3, PET_data_4)We read the PRE (Monthly Precipitation) data and prepare it.
# Open the first netCDF file with monthly Precipitation dataand store it
PRE <- nc_open("Data/cru_ts4.03.1981.1990.pre.dat.nc")
PRE_data <- ncvar_get(PRE, attributes(PRE$var)$names[1])
# Subset the data for the relevant time steps
PRE_data<-PRE_data[,,7:120]
# Close the file
nc_close(PRE)
# Open the second netCDF file with monthly Precipitation data and store it
PRE<-nc_open("Data/cru_ts4.03.1991.2000.pre.dat.nc")
PRE_data_2 <- ncvar_get(PRE, attributes(PRE$var)$names[1])
# Close the file
nc_close(PRE)
# Open the third netCDF file with monthly Precipitation data and store it
PRE<-nc_open("Data/cru_ts4.03.2001.2010.pre.dat.nc")
PRE_data_3 <- ncvar_get(PRE, attributes(PRE$var)$names[1])
# Close the file
nc_close(PRE)
# Open the fourth netCDF file with monthly Precipitation data and store it
PRE<-nc_open("Data/cru_ts4.03.2011.2018.pre.dat.nc")
PRE_data_4 <- ncvar_get(PRE, attributes(PRE$var)$names[1])
# Subset the data for the relevant time steps
PRE_data_4<-PRE_data_4[,,1:48]
# Get the longitude and latitude values
PRE_lon <- ncvar_get(PRE, attributes(PRE$dim)$names[1])
PRE_lat <- ncvar_get(PRE, attributes(PRE$dim)$names[2])
# Print a message to check if the longitudes and latitudes match with the ones from SPEI
print(paste(
"Check on longitudes",
all(PRE_lon == SPEI_lon),
"check on latitudes",
all(PRE_lat == SPEI_lat)
))## [1] "Check on longitudes TRUE check on latitudes TRUE"# Close the file
nc_close(PRE)
# Merge the 4 Precipitation datasets together
PRE<-array(0, dim=c(720,360,402))
PRE[,,]<-c(PRE_data, PRE_data_2,PRE_data_3,PRE_data_4)We clean-up the workspace by removing unnecessary objects.
# Remove the variables SPEI and SMoist from the workspace
rm(SPEI, SMoist)
# Rename the variables SPEI and Smoist variables
SPEI <- SPEI_relevant_years
SMoist <- SMoist_relevant_years
# Assign the variables lats and lons to the latitude and longitude of SPEI
lats <- SPEI_lat
lons <- SPEI_lon
# Remove variables that won't be used again from the workspace
rm(
SPEI_relevant_years,
SMoist_relevant_years,
SPEI_lat,
SMoist_lat,
SPEI_lon,
SMoist_lon,
SMoist_ordered,
NDVI_lon,
NDVI_lat,
SPEI_data,
SMoist_data,
NDVI_read,
PET_lat,
PET_lon,
PRE_lat,
PRE_lon,
PET_data,
PET_data_2,
PET_data_3,
PET_data_4,
PRE_data,
PRE_data_2,
PRE_data_3,
PRE_data_4
)We create arrays containing the latitudes and longitudes of the region (Continental Italy) we’ll analyze.
reg_latitudes <-
c(
46.75,
46.25,
46.25,
45.75,
45.25,
44.75,
44.25,
43.75,
43.75,
43.25,
42.75,
42.25,
41.75,
41.25,
40.75,
40.25,
39.75,
39.25,
38.75,
38.25
)
reg_long_min <-
c(
10.25,
8.25,
9.25,
6.75,
6.75,
6.75,
6.75,
7.25,
10.25,
10.25,
10.75,
11.75,
12.25,
13.25,
14.25,
15.25,
15.75,
16.25,
16.25,
15.75
)
reg_long_max <-
c(
12.75,
8.25,
13.75,
13.75,
12.25,
12.25,
12.25,
7.75,
13.75,
13.75,
14.25,
14.25,
15.75,
16.75,
17.75,
18.25,
16.25,
16.75,
16.75,
16.25
)We now scale the data and then format it as needed.
test_df <- create_test_df(reg_latitudes,
reg_long_min,
reg_long_max,
lats,
lons,
SPEI,
SMoist,
NDVI)
#We store the scaling factors for future re-use
scale_SMoist <- sd(test_df$SMoist)
scale_NDVI <- sd(test_df$NDVI, na.rm = TRUE)
#Since we don't scale the SPEI value we set the scale factor to 1
scale_SPEI <- 1
#We scale the data
test_df$SMoist <- (test_df$SMoist) / (scale_SMoist)
test_df$NDVI <- (test_df$NDVI) / (scale_NDVI)
#We create a dataframe for each of the variables
test_df_spei <-
data.frame(split(test_df$SPEI, as.factor(test_df$index)))
test_df_ndvi <-
data.frame(split(test_df$NDVI, as.factor(test_df$index)))
test_df_smoist <-
data.frame(split(test_df$SMoist, as.factor(test_df$index)))We create two other objects needed in the future analysis (the regions indexes and the weights matrix).
#index_pairs contains the Longitudes and Latitudes indexes of all the regions under consideration
index_pairs <- unique(test_df$index)
#weights is the connectivity matrix
weights <- prepare_weights(index_pairs, lats, lons)Unit Root testing
We now perform unit root testing as described in Section 3.3. We first look at the estimated \(\lambda\) for the SPEI variable on our data.
compute_lambda(test_df_spei, weights)## [1] 0.9071363We then look at the estimated \(\lambda\) for the SMoist variable.
compute_lambda(test_df_smoist, weights)## [1] 0.8726805We also look at the estimated \(\lambda\) for the NDVI variable.
compute_lambda(test_df_ndvi, weights)## [1] 0.3985638Finally, we perform a Monte Carlo simulation to obtain the critical values of our test.
lambda_mc <-
mc_critical_values(10000, 139, 500, weights, 3061792, FALSE)## [1] "We reject in the test 1 if lambda < 0.991387964490184"
## [1] "We reject in the test 2 if lambda < 0.996196614372254"
## [1] "We reject in the test 3 if lambda < 1.02291769144339"
## [1] "We reject in the test 4 if lambda < 1.05695372345613"
## [1] "We reject in the test 5 if lambda < 1.09510614689225"
## [1] "We reject in the test 6 if lambda < 1.13546696878769"
## [1] "We reject in the test 7 if lambda < 1.17463616654881"
## [1] "We reject in the test 8 if lambda < 1.20721156481164"
## [1] "We reject in the test 9 if lambda < 1.2234286905795"
## [1] "We reject in the test 10 if lambda < 1.2067064741227"Model fitting
We fit our model using SUR (implemented in the systemfit package) after having prepared the data for it.
spvar_df <-
prepare_spvar_df(test_df, test_df_spei, test_df_smoist, test_df_ndvi, lags =
1)
spvar_fit <-
systemfit(
c(
eqspei = spei ~ 0 + . -spei -ndvi -smoist,
eqndvi = ndvi ~ 0 + . -spei -ndvi -smoist,
eqsmoist = smoist ~ 0 + . -spei -ndvi -smoist
),
data = spvar_df,
method = "SUR"
)We also fit a model with two temporal lags to compare them.
spvar_df_2 <-
prepare_spvar_df(test_df, test_df_spei, test_df_smoist, test_df_ndvi, lags =
2)
spvar_fit_2 <-
systemfit(
c(
eqspei = spei ~ 0 + . -spei -ndvi -smoist,
eqndvi = ndvi ~ 0 + . -spei -ndvi -smoist,
eqsmoist = smoist ~ 0 + . -spei -ndvi -smoist
),
data = spvar_df_2,
method = "SUR"
)To make this comparison and choose which one to use we use BIC.
print(paste( "The BIC of the first model is", BIC(spvar_fit)))## [1] "The BIC of the first model is 131207.193675043"print(paste( "The BIC of the second model is", BIC(spvar_fit_2)))## [1] "The BIC of the second model is 143753.030804394"Since the BIC of the first model is substantially lower we choose it. We now look at the coefficients of this chosen model.
summary(spvar_fit)##
## systemfit results
## method: SUR
##
## N DF SSR detRCov OLS-R2 McElroy-R2
## system 145287 144852 53020.1 0.010642 0.617155 NA
##
## N DF SSR MSE RMSE R2 Adj R2
## eqspei 48789 48644 6872.68 0.141285 0.375879 0.840740 0.840269
## eqndvi 47709 47564 41672.47 0.876135 0.936021 0.120492 0.117830
## eqsmoist 48789 48644 4474.95 0.091994 0.303305 0.906683 0.906407
##
## The covariance matrix of the residuals used for estimation
## eqspei eqndvi eqsmoist
## eqspei 0.14142911 -0.0227940 0.00510905
## eqndvi -0.02279404 0.8761230 -0.07153691
## eqsmoist 0.00510905 -0.0715369 0.09216100
##
## The covariance matrix of the residuals
## eqspei eqndvi eqsmoist
## eqspei 0.14142885 -0.0227889 0.00510902
## eqndvi -0.02278892 0.8761347 -0.07155297
## eqsmoist 0.00510902 -0.0715530 0.09216346
##
## The correlations of the residuals
## eqspei eqndvi eqsmoist
## eqspei 1.0000000 -0.0647382 0.044745
## eqndvi -0.0647382 1.0000000 -0.251801
## eqsmoist 0.0447450 -0.2518006 1.000000
##
##
## SUR estimates for 'eqspei' (equation 1)
## Model Formula: spei ~ 0 + (ndvi + smoist + fe.1 + fe.2 + fe.3 + fe.4 + fe.5 +
## fe.6 + fe.7 + fe.8 + fe.9 + fe.10 + fe.11 + fe.12 + fe.13 +
## fe.14 + fe.15 + fe.16 + fe.17 + fe.18 + fe.19 + fe.20 + fe.21 +
## fe.22 + fe.23 + fe.24 + fe.25 + fe.26 + fe.27 + fe.28 + fe.29 +
## fe.30 + fe.31 + fe.32 + fe.33 + fe.34 + fe.35 + fe.36 + fe.37 +
## fe.38 + fe.39 + fe.40 + fe.41 + fe.42 + fe.43 + fe.44 + fe.45 +
## fe.46 + fe.47 + fe.48 + fe.49 + fe.50 + fe.51 + fe.52 + fe.53 +
## fe.54 + fe.55 + fe.56 + fe.57 + fe.58 + fe.59 + fe.60 + fe.61 +
## fe.62 + fe.63 + fe.64 + fe.65 + fe.66 + fe.67 + fe.68 + fe.69 +
## fe.70 + fe.71 + fe.72 + fe.73 + fe.74 + fe.75 + fe.76 + fe.77 +
## fe.78 + fe.79 + fe.80 + fe.81 + fe.82 + fe.83 + fe.84 + fe.85 +
## fe.86 + fe.87 + fe.88 + fe.89 + fe.90 + fe.91 + fe.92 + fe.93 +
## fe.94 + fe.95 + fe.96 + fe.97 + fe.98 + fe.99 + fe.100 +
## fe.101 + fe.102 + fe.103 + fe.104 + fe.105 + fe.106 + fe.107 +
## fe.108 + fe.109 + fe.110 + fe.111 + fe.112 + fe.113 + fe.114 +
## fe.115 + fe.116 + fe.117 + fe.118 + fe.119 + fe.120 + fe.121 +
## fe.122 + fe.123 + fe.124 + fe.125 + fe.126 + fe.127 + fe.128 +
## fe.129 + fe.130 + fe.131 + fe.132 + fe.133 + fe.134 + fe.135 +
## fe.136 + fe.137 + fe.138 + fe.139 + spei_lag + ndvi_lag +
## smoist_lag + sp_spei + sp_ndvi + sp_smoist) - spei - ndvi -
## smoist
##
## Estimate Std. Error t value Pr(>|t|)
## fe.1 -0.01183545 0.02259684 -0.52377 0.60044381
## fe.2 -0.02262707 0.02203690 -1.02678 0.30452861
## fe.3 -0.03930379 0.02185573 -1.79833 0.07213104 .
## fe.4 -0.03464724 0.02189803 -1.58221 0.11360852
## fe.5 -0.04546455 0.02174939 -2.09038 0.03658863 *
## fe.6 -0.04069721 0.02159535 -1.88454 0.05949851 .
## fe.7 -0.04652168 0.02156243 -2.15753 0.03096891 *
## fe.8 -0.05870931 0.02161661 -2.71594 0.00661121 **
## fe.9 -0.06149561 0.02141683 -2.87137 0.00408876 **
## fe.10 -0.07710600 0.02200197 -3.50450 0.00045787 ***
## fe.11 -0.07749231 0.02205000 -3.51439 0.00044117 ***
## fe.12 -0.07795694 0.02149009 -3.62758 0.00028639 ***
## fe.13 -0.06723482 0.02095264 -3.20890 0.00133332 **
## fe.14 -0.05821449 0.02083883 -2.79356 0.00521520 **
## fe.15 -0.03351792 0.02153550 -1.55640 0.11961886
## fe.16 -0.06305272 0.02155620 -2.92504 0.00344572 **
## fe.17 -0.07262124 0.02179201 -3.33247 0.00086142 ***
## fe.18 -0.09100797 0.02222406 -4.09502 4.2281e-05 ***
## fe.19 -0.09544689 0.02241314 -4.25852 2.0617e-05 ***
## fe.20 -0.09263960 0.02214717 -4.18291 2.8831e-05 ***
## fe.21 -0.08390869 0.02160800 -3.88322 0.00010322 ***
## fe.22 -0.06571049 0.02093290 -3.13910 0.00169568 **
## fe.23 -0.05389713 0.02097544 -2.56954 0.01018646 *
## fe.24 -0.06313391 0.02156410 -2.92773 0.00341602 **
## fe.25 -0.07319709 0.02202954 -3.32268 0.00089223 ***
## fe.26 -0.08625032 0.02260214 -3.81602 0.00013579 ***
## fe.27 -0.09131502 0.02237140 -4.08178 4.4764e-05 ***
## fe.28 -0.10294516 0.02254870 -4.56546 4.9962e-06 ***
## fe.29 -0.10172229 0.02231996 -4.55746 5.1902e-06 ***
## fe.30 -0.08339154 0.02134283 -3.90724 9.3483e-05 ***
## fe.31 -0.06610976 0.02088618 -3.16524 0.00155050 **
## fe.32 -0.06163741 0.02158465 -2.85561 0.00429720 **
## fe.33 -0.07392344 0.02235515 -3.30677 0.00094445 ***
## fe.34 -0.08405280 0.02306361 -3.64439 0.00026830 ***
## fe.35 -0.09037754 0.02342383 -3.85836 0.00011430 ***
## fe.36 -0.09274141 0.02304804 -4.02383 5.7345e-05 ***
## fe.37 -0.08913074 0.02201931 -4.04784 5.1771e-05 ***
## fe.38 -0.09040799 0.02164445 -4.17696 2.9595e-05 ***
## fe.39 -0.07283610 0.02101528 -3.46586 0.00052899 ***
## fe.40 -0.05827551 0.02182718 -2.66986 0.00759082 **
## fe.41 -0.08432885 0.02333687 -3.61355 0.00030234 ***
## fe.42 -0.09149210 0.02385832 -3.83481 0.00012582 ***
## fe.43 -0.09561791 0.02376359 -4.02372 5.7374e-05 ***
## fe.44 -0.09580446 0.02370435 -4.04164 5.3160e-05 ***
## fe.45 -0.08525938 0.02264917 -3.76435 0.00016718 ***
## fe.46 -0.06628662 0.02179106 -3.04192 0.00235201 **
## fe.47 -0.07446181 0.02274654 -3.27354 0.00106282 **
## fe.48 -0.08637753 0.02400696 -3.59802 0.00032097 ***
## fe.49 -0.09219954 0.02464140 -3.74165 0.00018303 ***
## fe.50 -0.09588738 0.02450102 -3.91361 9.1051e-05 ***
## fe.51 -0.09946895 0.02413193 -4.12188 3.7641e-05 ***
## fe.52 -0.09599468 0.02400490 -3.99896 6.3715e-05 ***
## fe.53 -0.05947229 0.02348294 -2.53257 0.01132597 *
## fe.54 -0.04477553 0.02360326 -1.89701 0.05783308 .
## fe.55 -0.08141773 0.02427065 -3.35458 0.00079548 ***
## fe.56 -0.08864370 0.02453237 -3.61334 0.00030259 ***
## fe.57 -0.09287429 0.02499849 -3.71520 0.00020327 ***
## fe.58 -0.09230880 0.02491885 -3.70438 0.00021215 ***
## fe.59 -0.08791247 0.02429194 -3.61900 0.00029605 ***
## fe.60 -0.08990450 0.02398859 -3.74780 0.00017860 ***
## fe.61 -0.08271892 0.02298823 -3.59832 0.00032060 ***
## fe.62 -0.06361303 0.02147161 -2.96266 0.00305143 **
## fe.63 -0.06611646 0.02137874 -3.09263 0.00198505 **
## fe.64 -0.07980833 0.02379633 -3.35381 0.00079768 ***
## fe.65 -0.08925187 0.02490797 -3.58327 0.00033966 ***
## fe.66 -0.09278291 0.02537060 -3.65710 0.00025535 ***
## fe.67 -0.09085287 0.02514619 -3.61299 0.00030300 ***
## fe.68 -0.09225742 0.02476734 -3.72496 0.00019556 ***
## fe.69 -0.07438808 0.02321059 -3.20492 0.00135187 **
## fe.70 -0.07956095 0.02298645 -3.46121 0.00053821 ***
## fe.71 -0.04698539 0.02371432 -1.98131 0.04756227 *
## fe.72 -0.05720536 0.02154458 -2.65521 0.00792852 **
## fe.73 -0.08914541 0.02248838 -3.96407 7.3789e-05 ***
## fe.74 -0.09311164 0.02299749 -4.04877 5.1566e-05 ***
## fe.75 -0.09047707 0.02284122 -3.96113 7.4702e-05 ***
## fe.76 -0.08452781 0.02306750 -3.66437 0.00024821 ***
## fe.77 -0.08322574 0.02342017 -3.55359 0.00038037 ***
## fe.78 -0.07847685 0.02373804 -3.30595 0.00094722 ***
## fe.79 -0.09232268 0.02532833 -3.64504 0.00026763 ***
## fe.80 -0.09253221 0.02560068 -3.61444 0.00030130 ***
## fe.81 -0.09503188 0.02600620 -3.65420 0.00025826 ***
## fe.82 -0.09138468 0.02488921 -3.67166 0.00024124 ***
## fe.83 -0.08505867 0.02396912 -3.54868 0.00038754 ***
## fe.84 -0.05931886 0.02200658 -2.69551 0.00703059 **
## fe.85 -0.08657876 0.02313011 -3.74312 0.00018196 ***
## fe.86 -0.09904682 0.02386416 -4.15044 3.3240e-05 ***
## fe.87 -0.10404560 0.02450323 -4.24620 2.1783e-05 ***
## fe.88 -0.09895034 0.02471956 -4.00292 6.2659e-05 ***
## fe.89 -0.10528932 0.02496481 -4.21751 2.4746e-05 ***
## fe.90 -0.10406835 0.02574255 -4.04266 5.2930e-05 ***
## fe.91 -0.10016060 0.02582251 -3.87881 0.00010511 ***
## fe.92 -0.10211773 0.02658243 -3.84155 0.00012241 ***
## fe.93 -0.10232607 0.02697612 -3.79321 0.00014889 ***
## fe.94 -0.09516635 0.02567061 -3.70721 0.00020979 ***
## fe.95 -0.08085568 0.02380964 -3.39592 0.00068452 ***
## fe.96 -0.05054149 0.02197040 -2.30044 0.02142774 *
## fe.97 -0.07664014 0.02336601 -3.27998 0.00103886 **
## fe.98 -0.08772090 0.02398581 -3.65720 0.00025525 ***
## fe.99 -0.09590126 0.02483334 -3.86179 0.00011270 ***
## fe.100 -0.09761283 0.02524094 -3.86724 0.00011022 ***
## fe.101 -0.10447717 0.02631025 -3.97097 7.1684e-05 ***
## fe.102 -0.10685948 0.02661054 -4.01568 5.9363e-05 ***
## fe.103 -0.10377297 0.02673175 -3.88201 0.00010373 ***
## fe.104 -0.09674882 0.02659154 -3.63833 0.00027469 ***
## fe.105 -0.09231261 0.02621846 -3.52090 0.00043048 ***
## fe.106 -0.08917392 0.02558100 -3.48594 0.00049084 ***
## fe.107 -0.07811449 0.02392870 -3.26447 0.00109745 **
## fe.108 -0.02082546 0.02286359 -0.91086 0.36237529
## fe.109 -0.04485275 0.02218745 -2.02154 0.04322970 *
## fe.110 -0.05909096 0.02263084 -2.61108 0.00902842 **
## fe.111 -0.07932768 0.02432752 -3.26082 0.00111167 **
## fe.112 -0.07144802 0.02350738 -3.03939 0.00237186 **
## fe.113 -0.08505697 0.02484408 -3.42363 0.00061841 ***
## fe.114 -0.09034513 0.02580058 -3.50167 0.00046276 ***
## fe.115 -0.09571612 0.02602492 -3.67786 0.00023545 ***
## fe.116 -0.09486270 0.02601958 -3.64582 0.00026682 ***
## fe.117 -0.08852947 0.02553935 -3.46640 0.00052794 ***
## fe.118 -0.09408643 0.02554357 -3.68337 0.00023042 ***
## fe.119 -0.09340345 0.02484120 -3.76002 0.00017010 ***
## fe.120 -0.08414831 0.02324697 -3.61975 0.00029518 ***
## fe.121 -0.06943038 0.02264973 -3.06540 0.00217502 **
## fe.122 -0.04809657 0.02259785 -2.12837 0.03331147 *
## fe.123 -0.01762394 0.02545126 -0.69246 0.48865277
## fe.124 -0.04590827 0.02354873 -1.94950 0.05124142 .
## fe.125 -0.06015963 0.02354141 -2.55548 0.01060711 *
## fe.126 -0.08391095 0.02526263 -3.32154 0.00089587 ***
## fe.127 -0.07356572 0.02441697 -3.01289 0.00258903 **
## fe.128 -0.07807990 0.02514235 -3.10551 0.00190058 **
## fe.129 -0.08406474 0.02463689 -3.41215 0.00064505 ***
## fe.130 -0.08639713 0.02445491 -3.53292 0.00041139 ***
## fe.131 -0.08685821 0.02380544 -3.64867 0.00026388 ***
## fe.132 -0.06610577 0.02289880 -2.88687 0.00389271 **
## fe.133 -0.04939431 0.02356808 -2.09581 0.03610384 *
## fe.134 -0.03457177 0.02218661 -1.55823 0.11918595
## fe.135 -0.05478925 0.02279275 -2.40380 0.01622927 *
## fe.136 -0.06035373 0.02287976 -2.63787 0.00834562 **
## fe.137 -0.06129471 0.02269938 -2.70028 0.00693047 **
## fe.138 -0.05798174 0.02283436 -2.53923 0.01111269 *
## fe.139 -0.04568695 0.02318442 -1.97059 0.04877662 *
## spei_lag 0.94191511 0.00548691 171.66579 < 2.22e-16 ***
## ndvi_lag -0.00344450 0.00395246 -0.87148 0.38349521
## smoist_lag -0.00739863 0.00515436 -1.43541 0.15117613
## sp_spei -0.04201634 0.00721678 -5.82204 5.8499e-09 ***
## sp_ndvi 0.01682692 0.00549154 3.06415 0.00218407 **
## sp_smoist 0.02404245 0.00759167 3.16695 0.00154141 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.375879 on 48644 degrees of freedom
## Number of observations: 48789 Degrees of Freedom: 48644
## SSR: 6872.676626 MSE: 0.141285 Root MSE: 0.375879
## Multiple R-Squared: 0.84074 Adjusted R-Squared: 0.840269
##
##
## SUR estimates for 'eqndvi' (equation 2)
## Model Formula: ndvi ~ 0 + (spei + smoist + fe.1 + fe.2 + fe.3 + fe.4 + fe.5 +
## fe.6 + fe.7 + fe.8 + fe.9 + fe.10 + fe.11 + fe.12 + fe.13 +
## fe.14 + fe.15 + fe.16 + fe.17 + fe.18 + fe.19 + fe.20 + fe.21 +
## fe.22 + fe.23 + fe.24 + fe.25 + fe.26 + fe.27 + fe.28 + fe.29 +
## fe.30 + fe.31 + fe.32 + fe.33 + fe.34 + fe.35 + fe.36 + fe.37 +
## fe.38 + fe.39 + fe.40 + fe.41 + fe.42 + fe.43 + fe.44 + fe.45 +
## fe.46 + fe.47 + fe.48 + fe.49 + fe.50 + fe.51 + fe.52 + fe.53 +
## fe.54 + fe.55 + fe.56 + fe.57 + fe.58 + fe.59 + fe.60 + fe.61 +
## fe.62 + fe.63 + fe.64 + fe.65 + fe.66 + fe.67 + fe.68 + fe.69 +
## fe.70 + fe.71 + fe.72 + fe.73 + fe.74 + fe.75 + fe.76 + fe.77 +
## fe.78 + fe.79 + fe.80 + fe.81 + fe.82 + fe.83 + fe.84 + fe.85 +
## fe.86 + fe.87 + fe.88 + fe.89 + fe.90 + fe.91 + fe.92 + fe.93 +
## fe.94 + fe.95 + fe.96 + fe.97 + fe.98 + fe.99 + fe.100 +
## fe.101 + fe.102 + fe.103 + fe.104 + fe.105 + fe.106 + fe.107 +
## fe.108 + fe.109 + fe.110 + fe.111 + fe.112 + fe.113 + fe.114 +
## fe.115 + fe.116 + fe.117 + fe.118 + fe.119 + fe.120 + fe.121 +
## fe.122 + fe.123 + fe.124 + fe.125 + fe.126 + fe.127 + fe.128 +
## fe.129 + fe.130 + fe.131 + fe.132 + fe.133 + fe.134 + fe.135 +
## fe.136 + fe.137 + fe.138 + fe.139 + spei_lag + ndvi_lag +
## smoist_lag + sp_spei + sp_ndvi + sp_smoist) - spei - ndvi -
## smoist
##
## Estimate Std. Error t value Pr(>|t|)
## fe.1 1.2120519 0.0569747 21.27349 < 2.22e-16 ***
## fe.2 1.0573564 0.0555474 19.03522 < 2.22e-16 ***
## fe.3 1.2663599 0.0550890 22.98753 < 2.22e-16 ***
## fe.4 0.8753603 0.0551325 15.87740 < 2.22e-16 ***
## fe.5 1.2647504 0.0547700 23.09202 < 2.22e-16 ***
## fe.6 1.1302447 0.0543893 20.78066 < 2.22e-16 ***
## fe.7 1.1581469 0.0542943 21.33091 < 2.22e-16 ***
## fe.8 1.0843162 0.0543759 19.94112 < 2.22e-16 ***
## fe.9 1.1854510 0.0539368 21.97852 < 2.22e-16 ***
## fe.10 1.0900740 0.0553557 19.69219 < 2.22e-16 ***
## fe.11 1.0196376 0.0554796 18.37860 < 2.22e-16 ***
## fe.12 0.9001672 0.0540747 16.64672 < 2.22e-16 ***
## fe.13 0.8833363 0.0527854 16.73447 < 2.22e-16 ***
## fe.14 0.9385153 0.0524293 17.90060 < 2.22e-16 ***
## fe.15 0.9467465 0.0541071 17.49763 < 2.22e-16 ***
## fe.16 0.9951563 0.0543015 18.32648 < 2.22e-16 ***
## fe.17 1.1828420 0.0548921 21.54849 < 2.22e-16 ***
## fe.18 1.1214414 0.0559187 20.05485 < 2.22e-16 ***
## fe.19 1.0828807 0.0563936 19.20218 < 2.22e-16 ***
## fe.20 0.9833455 0.0557243 17.64661 < 2.22e-16 ***
## fe.21 1.0206563 0.0543670 18.77345 < 2.22e-16 ***
## fe.22 0.9985327 0.0526650 18.96008 < 2.22e-16 ***
## fe.23 1.0322180 0.0526981 19.58739 < 2.22e-16 ***
## fe.24 1.1546038 0.0542671 21.27630 < 2.22e-16 ***
## fe.25 1.2572113 0.0555047 22.65054 < 2.22e-16 ***
## fe.26 1.2546151 0.0569460 22.03168 < 2.22e-16 ***
## fe.27 1.1008521 0.0563599 19.53256 < 2.22e-16 ***
## fe.28 0.9786703 0.0567391 17.24860 < 2.22e-16 ***
## fe.29 0.9759688 0.0561592 17.37861 < 2.22e-16 ***
## fe.30 0.9275800 0.0536996 17.27350 < 2.22e-16 ***
## fe.31 0.8636603 0.0525443 16.43680 < 2.22e-16 ***
## fe.32 1.1435120 0.0542670 21.07195 < 2.22e-16 ***
## fe.33 1.2018609 0.0562747 21.35704 < 2.22e-16 ***
## fe.34 1.1633356 0.0580668 20.03445 < 2.22e-16 ***
## fe.35 1.0220237 0.0590290 17.31392 < 2.22e-16 ***
## fe.36 1.1125924 0.0580701 19.15947 < 2.22e-16 ***
## fe.37 1.1057843 0.0554102 19.95634 < 2.22e-16 ***
## fe.38 1.0071070 0.0544618 18.49200 < 2.22e-16 ***
## fe.39 1.0735827 0.0528694 20.30633 < 2.22e-16 ***
## fe.40 1.2768473 0.0548771 23.26740 < 2.22e-16 ***
## fe.41 1.3264694 0.0587044 22.59572 < 2.22e-16 ***
## fe.42 1.2875222 0.0600708 21.43341 < 2.22e-16 ***
## fe.43 1.0018790 0.0598354 16.74393 < 2.22e-16 ***
## fe.44 1.0063588 0.0597460 16.84396 < 2.22e-16 ***
## fe.45 1.1099621 0.0570714 19.44866 < 2.22e-16 ***
## fe.46 1.2467130 0.0547796 22.75870 < 2.22e-16 ***
## fe.47 1.2427840 0.0572030 21.72585 < 2.22e-16 ***
## fe.48 1.2051098 0.0603996 19.95228 < 2.22e-16 ***
## fe.49 1.2377734 0.0620064 19.96204 < 2.22e-16 ***
## fe.50 1.2177184 0.0616932 19.73829 < 2.22e-16 ***
## fe.51 1.0027907 0.0607534 16.50591 < 2.22e-16 ***
## fe.52 1.0479241 0.0604876 17.32462 < 2.22e-16 ***
## fe.53 0.6359064 0.0601081 10.57938 < 2.22e-16 ***
## fe.54 1.0492095 0.0596337 17.59425 < 2.22e-16 ***
## fe.55 1.3473167 0.0610542 22.06755 < 2.22e-16 ***
## fe.56 1.1930170 0.0617281 19.32697 < 2.22e-16 ***
## fe.57 1.0879523 0.0629063 17.29480 < 2.22e-16 ***
## fe.58 1.1075361 0.0627078 17.66184 < 2.22e-16 ***
## fe.59 1.0867516 0.0611688 17.76644 < 2.22e-16 ***
## fe.60 1.0509511 0.0604556 17.38385 < 2.22e-16 ***
## fe.61 0.9150440 0.0579258 15.79684 < 2.22e-16 ***
## fe.62 1.1485239 0.0540385 21.25380 < 2.22e-16 ***
## fe.63 1.1310983 0.0537337 21.05007 < 2.22e-16 ***
## fe.64 1.3061491 0.0599063 21.80320 < 2.22e-16 ***
## fe.65 1.2382424 0.0626647 19.75982 < 2.22e-16 ***
## fe.66 1.2052923 0.0638440 18.87871 < 2.22e-16 ***
## fe.67 1.1632090 0.0632715 18.38441 < 2.22e-16 ***
## fe.68 1.1180571 0.0623166 17.94157 < 2.22e-16 ***
## fe.69 1.0072137 0.0584318 17.23744 < 2.22e-16 ***
## fe.70 0.8800142 0.0579269 15.19179 < 2.22e-16 ***
## fe.71 0.6386129 0.0604971 10.55609 < 2.22e-16 ***
## fe.72 0.6834496 0.0541428 12.62310 < 2.22e-16 ***
## fe.73 0.8518486 0.0565441 15.06521 < 2.22e-16 ***
## fe.74 1.0247893 0.0578849 17.70390 < 2.22e-16 ***
## fe.75 1.1326726 0.0575054 19.69680 < 2.22e-16 ***
## fe.76 0.9910942 0.0580197 17.08203 < 2.22e-16 ***
## fe.77 1.2579619 0.0589011 21.35718 < 2.22e-16 ***
## fe.78 1.3052090 0.0597061 21.86056 < 2.22e-16 ***
## fe.79 1.1630878 0.0637298 18.25031 < 2.22e-16 ***
## fe.80 1.1507491 0.0644418 17.85718 < 2.22e-16 ***
## fe.81 1.2434039 0.0654427 18.99987 < 2.22e-16 ***
## fe.82 1.1855180 0.0626236 18.93087 < 2.22e-16 ***
## fe.83 1.0034460 0.0602859 16.64478 < 2.22e-16 ***
## fe.84 0.5967273 0.0553697 10.77714 < 2.22e-16 ***
## fe.85 1.0108272 0.0581612 17.37976 < 2.22e-16 ***
## fe.86 1.2469102 0.0600156 20.77644 < 2.22e-16 ***
## fe.87 1.2136059 0.0616283 19.69233 < 2.22e-16 ***
## fe.88 1.1672299 0.0621848 18.77035 < 2.22e-16 ***
## fe.89 1.1497458 0.0628152 18.30363 < 2.22e-16 ***
## fe.90 1.1909007 0.0647692 18.38683 < 2.22e-16 ***
## fe.91 1.1937804 0.0649760 18.37264 < 2.22e-16 ***
## fe.92 1.1795183 0.0668945 17.63251 < 2.22e-16 ***
## fe.93 1.0756458 0.0678874 15.84456 < 2.22e-16 ***
## fe.94 1.0858226 0.0645903 16.81093 < 2.22e-16 ***
## fe.95 0.8448379 0.0598904 14.10639 < 2.22e-16 ***
## fe.96 0.5397137 0.0552077 9.77606 < 2.22e-16 ***
## fe.97 0.7938626 0.0587495 13.51266 < 2.22e-16 ***
## fe.98 1.2260813 0.0603080 20.33033 < 2.22e-16 ***
## fe.99 1.1430462 0.0624578 18.30110 < 2.22e-16 ***
## fe.100 1.1593143 0.0634965 18.25792 < 2.22e-16 ***
## fe.101 1.2201937 0.0661989 18.43223 < 2.22e-16 ***
## fe.102 1.2606224 0.0669638 18.82544 < 2.22e-16 ***
## fe.103 1.2750852 0.0672682 18.95524 < 2.22e-16 ***
## fe.104 1.2396389 0.0669057 18.52815 < 2.22e-16 ***
## fe.105 1.1556115 0.0659681 17.51772 < 2.22e-16 ***
## fe.106 1.1217694 0.0642955 17.44710 < 2.22e-16 ***
## fe.107 0.7448362 0.0601211 12.38894 < 2.22e-16 ***
## fe.108 0.6079940 0.0574346 10.58584 < 2.22e-16 ***
## fe.109 0.5766444 0.0557592 10.34168 < 2.22e-16 ***
## fe.110 0.6695410 0.0568977 11.76744 < 2.22e-16 ***
## fe.111 1.1358543 0.0611818 18.56524 < 2.22e-16 ***
## fe.112 1.1615647 0.0591084 19.65142 < 2.22e-16 ***
## fe.113 1.1005487 0.0624916 17.61116 < 2.22e-16 ***
## fe.114 1.0925548 0.0649083 16.83228 < 2.22e-16 ***
## fe.115 1.0111831 0.0654748 15.44386 < 2.22e-16 ***
## fe.116 1.0604944 0.0654582 16.20110 < 2.22e-16 ***
## fe.117 1.0519830 0.0642522 16.37271 < 2.22e-16 ***
## fe.118 1.1992921 0.0642442 18.66772 < 2.22e-16 ***
## fe.119 1.2040892 0.0624212 19.28975 < 2.22e-16 ***
## fe.120 1.0916963 0.0584469 18.67842 < 2.22e-16 ***
## fe.121 1.0283254 0.0569921 18.04331 < 2.22e-16 ***
## fe.122 1.1576583 0.0567850 20.38670 < 2.22e-16 ***
## fe.123 0.5785537 0.0639908 9.04121 < 2.22e-16 ***
## fe.124 0.7564577 0.0591920 12.77973 < 2.22e-16 ***
## fe.125 0.5793884 0.0591915 9.78837 < 2.22e-16 ***
## fe.126 0.6604226 0.0635289 10.39563 < 2.22e-16 ***
## fe.127 0.7685322 0.0614054 12.51572 < 2.22e-16 ***
## fe.128 1.1888408 0.0631734 18.81869 < 2.22e-16 ***
## fe.129 0.7805692 0.0619057 12.60901 < 2.22e-16 ***
## fe.130 0.9217613 0.0614434 15.00179 < 2.22e-16 ***
## fe.131 1.0431865 0.0598535 17.42899 < 2.22e-16 ***
## fe.132 1.0195162 0.0575432 17.71739 < 2.22e-16 ***
## fe.133 0.8589186 0.0592140 14.50534 < 2.22e-16 ***
## fe.134 0.4172274 0.0557446 7.48462 7.2831e-14 ***
## fe.135 0.4697165 0.0572543 8.20404 2.2204e-16 ***
## fe.136 0.8320678 0.0574599 14.48085 < 2.22e-16 ***
## fe.137 0.7570454 0.0570026 13.28089 < 2.22e-16 ***
## fe.138 0.5982570 0.0573352 10.43438 < 2.22e-16 ***
## fe.139 0.6151329 0.0581944 10.57031 < 2.22e-16 ***
## spei_lag -0.0100348 0.0137880 -0.72779 0.4667437
## ndvi_lag 0.0416286 0.0099459 4.18550 2.8504e-05 ***
## smoist_lag 0.0347218 0.0129525 2.68071 0.0073491 **
## sp_spei 0.0588650 0.0181155 3.24944 0.0011571 **
## sp_ndvi 0.2895066 0.0138422 20.91478 < 2.22e-16 ***
## sp_smoist -0.1112368 0.0190841 -5.82878 5.6195e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.936021 on 47564 degrees of freedom
## Number of observations: 47709 Degrees of Freedom: 47564
## SSR: 41672.469808 MSE: 0.876135 Root MSE: 0.936021
## Multiple R-Squared: 0.120492 Adjusted R-Squared: 0.11783
##
##
## SUR estimates for 'eqsmoist' (equation 3)
## Model Formula: smoist ~ 0 + (spei + ndvi + fe.1 + fe.2 + fe.3 + fe.4 + fe.5 +
## fe.6 + fe.7 + fe.8 + fe.9 + fe.10 + fe.11 + fe.12 + fe.13 +
## fe.14 + fe.15 + fe.16 + fe.17 + fe.18 + fe.19 + fe.20 + fe.21 +
## fe.22 + fe.23 + fe.24 + fe.25 + fe.26 + fe.27 + fe.28 + fe.29 +
## fe.30 + fe.31 + fe.32 + fe.33 + fe.34 + fe.35 + fe.36 + fe.37 +
## fe.38 + fe.39 + fe.40 + fe.41 + fe.42 + fe.43 + fe.44 + fe.45 +
## fe.46 + fe.47 + fe.48 + fe.49 + fe.50 + fe.51 + fe.52 + fe.53 +
## fe.54 + fe.55 + fe.56 + fe.57 + fe.58 + fe.59 + fe.60 + fe.61 +
## fe.62 + fe.63 + fe.64 + fe.65 + fe.66 + fe.67 + fe.68 + fe.69 +
## fe.70 + fe.71 + fe.72 + fe.73 + fe.74 + fe.75 + fe.76 + fe.77 +
## fe.78 + fe.79 + fe.80 + fe.81 + fe.82 + fe.83 + fe.84 + fe.85 +
## fe.86 + fe.87 + fe.88 + fe.89 + fe.90 + fe.91 + fe.92 + fe.93 +
## fe.94 + fe.95 + fe.96 + fe.97 + fe.98 + fe.99 + fe.100 +
## fe.101 + fe.102 + fe.103 + fe.104 + fe.105 + fe.106 + fe.107 +
## fe.108 + fe.109 + fe.110 + fe.111 + fe.112 + fe.113 + fe.114 +
## fe.115 + fe.116 + fe.117 + fe.118 + fe.119 + fe.120 + fe.121 +
## fe.122 + fe.123 + fe.124 + fe.125 + fe.126 + fe.127 + fe.128 +
## fe.129 + fe.130 + fe.131 + fe.132 + fe.133 + fe.134 + fe.135 +
## fe.136 + fe.137 + fe.138 + fe.139 + spei_lag + ndvi_lag +
## smoist_lag + sp_spei + sp_ndvi + sp_smoist) - spei - ndvi -
## smoist
##
## Estimate Std. Error t value Pr(>|t|)
## fe.1 0.25570974 0.01822538 14.03042 < 2e-16 ***
## fe.2 0.32663423 0.01777409 18.37699 < 2e-16 ***
## fe.3 0.41264187 0.01762800 23.40832 < 2e-16 ***
## fe.4 0.33970642 0.01766345 19.23217 < 2e-16 ***
## fe.5 0.43800145 0.01754331 24.96687 < 2e-16 ***
## fe.6 0.35724484 0.01741890 20.50904 < 2e-16 ***
## fe.7 0.39475841 0.01739260 22.69692 < 2e-16 ***
## fe.8 0.45192146 0.01743747 25.91669 < 2e-16 ***
## fe.9 0.40682996 0.01727497 23.55026 < 2e-16 ***
## fe.10 0.50900518 0.01774811 28.67940 < 2e-16 ***
## fe.11 0.49774811 0.01778679 27.98414 < 2e-16 ***
## fe.12 0.39350469 0.01733505 22.69995 < 2e-16 ***
## fe.13 0.33386160 0.01690017 19.75492 < 2e-16 ***
## fe.14 0.30394696 0.01680984 18.08149 < 2e-16 ***
## fe.15 0.21282476 0.01737341 12.25003 < 2e-16 ***
## fe.16 0.43158937 0.01738709 24.82240 < 2e-16 ***
## fe.17 0.48199884 0.01757737 27.42156 < 2e-16 ***
## fe.18 0.52091658 0.01792717 29.05738 < 2e-16 ***
## fe.19 0.51832788 0.01807971 28.66903 < 2e-16 ***
## fe.20 0.48548318 0.01786517 27.17485 < 2e-16 ***
## fe.21 0.41843741 0.01743026 24.00638 < 2e-16 ***
## fe.22 0.34084683 0.01688574 20.18548 < 2e-16 ***
## fe.23 0.27821955 0.01692163 16.44165 < 2e-16 ***
## fe.24 0.44648664 0.01739462 25.66809 < 2e-16 ***
## fe.25 0.51160409 0.01776866 28.79249 < 2e-16 ***
## fe.26 0.56989207 0.01823054 31.26029 < 2e-16 ***
## fe.27 0.54796164 0.01804453 30.36719 < 2e-16 ***
## fe.28 0.50844539 0.01818897 27.95349 < 2e-16 ***
## fe.29 0.47211562 0.01800455 26.22202 < 2e-16 ***
## fe.30 0.40480183 0.01721635 23.51264 < 2e-16 ***
## fe.31 0.32945217 0.01684812 19.55424 < 2e-16 ***
## fe.32 0.45367749 0.01741230 26.05501 < 2e-16 ***
## fe.33 0.53914532 0.01803237 29.89876 < 2e-16 ***
## fe.34 0.61685566 0.01860365 33.15778 < 2e-16 ***
## fe.35 0.61001954 0.01889305 32.28804 < 2e-16 ***
## fe.36 0.58260917 0.01859019 31.33961 < 2e-16 ***
## fe.37 0.48196875 0.01776186 27.13503 < 2e-16 ***
## fe.38 0.41838241 0.01745959 23.96289 < 2e-16 ***
## fe.39 0.35096711 0.01695226 20.70327 < 2e-16 ***
## fe.40 0.48363926 0.01760795 27.46710 < 2e-16 ***
## fe.41 0.61546170 0.01882514 32.69360 < 2e-16 ***
## fe.42 0.66506816 0.01924462 34.55865 < 2e-16 ***
## fe.43 0.68503948 0.01916814 35.73844 < 2e-16 ***
## fe.44 0.63242351 0.01911910 33.07811 < 2e-16 ***
## fe.45 0.53825590 0.01826834 29.46387 < 2e-16 ***
## fe.46 0.46635822 0.01757895 26.52936 < 2e-16 ***
## fe.47 0.54524317 0.01834929 29.71467 < 2e-16 ***
## fe.48 0.63657145 0.01936550 32.87142 < 2e-16 ***
## fe.49 0.68501725 0.01987706 34.46271 < 2e-16 ***
## fe.50 0.70747043 0.01976296 35.79780 < 2e-16 ***
## fe.51 0.69913353 0.01946546 35.91663 < 2e-16 ***
## fe.52 0.64114423 0.01936185 33.11378 < 2e-16 ***
## fe.53 0.53104549 0.01892196 28.06504 < 2e-16 ***
## fe.54 0.57930864 0.01903458 30.43454 < 2e-16 ***
## fe.55 0.65465172 0.01957839 33.43746 < 2e-16 ***
## fe.56 0.67455445 0.01978918 34.08703 < 2e-16 ***
## fe.57 0.70779311 0.02016507 35.09995 < 2e-16 ***
## fe.58 0.71857203 0.02010079 35.74844 < 2e-16 ***
## fe.59 0.70314751 0.01959427 35.88537 < 2e-16 ***
## fe.60 0.65471782 0.01934850 33.83817 < 2e-16 ***
## fe.61 0.56292224 0.01854182 30.35960 < 2e-16 ***
## fe.62 0.37616205 0.01731992 21.71847 < 2e-16 ***
## fe.63 0.39950996 0.01724651 23.16468 < 2e-16 ***
## fe.64 0.65346850 0.01919481 34.04403 < 2e-16 ***
## fe.65 0.71720526 0.02009235 35.69544 < 2e-16 ***
## fe.66 0.72429619 0.02046521 35.39158 < 2e-16 ***
## fe.67 0.73314297 0.02028436 36.14326 < 2e-16 ***
## fe.68 0.70863986 0.01997878 35.46962 < 2e-16 ***
## fe.69 0.63228993 0.01872232 33.77199 < 2e-16 ***
## fe.70 0.58492264 0.01854027 31.54877 < 2e-16 ***
## fe.71 0.52949611 0.01911237 27.70436 < 2e-16 ***
## fe.72 0.40940366 0.01738047 23.55539 < 2e-16 ***
## fe.73 0.54748979 0.01814124 30.17929 < 2e-16 ***
## fe.74 0.59248634 0.01855065 31.93884 < 2e-16 ***
## fe.75 0.58039704 0.01842430 31.50172 < 2e-16 ***
## fe.76 0.62440212 0.01860800 33.55557 < 2e-16 ***
## fe.77 0.67114704 0.01889262 35.52430 < 2e-16 ***
## fe.78 0.70983576 0.01914892 37.06923 < 2e-16 ***
## fe.79 0.78628394 0.02043129 38.48431 < 2e-16 ***
## fe.80 0.76779729 0.02065041 37.18073 < 2e-16 ***
## fe.81 0.72835511 0.02097794 34.72006 < 2e-16 ***
## fe.82 0.70321674 0.02007709 35.02583 < 2e-16 ***
## fe.83 0.60989259 0.01933536 31.54286 < 2e-16 ***
## fe.84 0.48860964 0.01775177 27.52455 < 2e-16 ***
## fe.85 0.59225085 0.01865885 31.74101 < 2e-16 ***
## fe.86 0.64749977 0.01925083 33.63491 < 2e-16 ***
## fe.87 0.68596365 0.01976624 34.70380 < 2e-16 ***
## fe.88 0.73879063 0.01994049 37.04978 < 2e-16 ***
## fe.89 0.78947154 0.02013804 39.20300 < 2e-16 ***
## fe.90 0.80617283 0.02076548 38.82274 < 2e-16 ***
## fe.91 0.81353815 0.02082986 39.05634 < 2e-16 ***
## fe.92 0.76362634 0.02144272 35.61239 < 2e-16 ***
## fe.93 0.70833361 0.02176024 32.55174 < 2e-16 ***
## fe.94 0.67043467 0.02070740 32.37658 < 2e-16 ***
## fe.95 0.57389141 0.01920659 29.87993 < 2e-16 ***
## fe.96 0.49001664 0.01772410 27.64691 < 2e-16 ***
## fe.97 0.58583901 0.01884925 31.08023 < 2e-16 ***
## fe.98 0.64303802 0.01934924 33.23324 < 2e-16 ***
## fe.99 0.70272544 0.02003255 35.07918 < 2e-16 ***
## fe.100 0.74608927 0.02036107 36.64294 < 2e-16 ***
## fe.101 0.78487513 0.02122339 36.98161 < 2e-16 ***
## fe.102 0.79232721 0.02146543 36.91179 < 2e-16 ***
## fe.103 0.79623657 0.02156321 36.92570 < 2e-16 ***
## fe.104 0.76740489 0.02145032 35.77592 < 2e-16 ***
## fe.105 0.72902179 0.02114934 34.47018 < 2e-16 ***
## fe.106 0.69076460 0.02063660 33.47280 < 2e-16 ***
## fe.107 0.60415465 0.01930411 31.29669 < 2e-16 ***
## fe.108 0.44639346 0.01844503 24.20129 < 2e-16 ***
## fe.109 0.53869297 0.01789908 30.09613 < 2e-16 ***
## fe.110 0.59947332 0.01825626 32.83659 < 2e-16 ***
## fe.111 0.70095682 0.01962459 35.71829 < 2e-16 ***
## fe.112 0.70596970 0.01896322 37.22837 < 2e-16 ***
## fe.113 0.79879983 0.02004107 39.85814 < 2e-16 ***
## fe.114 0.79708154 0.02081243 38.29834 < 2e-16 ***
## fe.115 0.79508778 0.02099336 37.87330 < 2e-16 ***
## fe.116 0.79469041 0.02098912 37.86202 < 2e-16 ***
## fe.117 0.77965837 0.02060168 37.84441 < 2e-16 ***
## fe.118 0.76114827 0.02060549 36.93911 < 2e-16 ***
## fe.119 0.72397267 0.02004009 36.12621 < 2e-16 ***
## fe.120 0.64746420 0.01875330 34.52535 < 2e-16 ***
## fe.121 0.60185764 0.01827050 32.94150 < 2e-16 ***
## fe.122 0.48913056 0.01823026 26.83069 < 2e-16 ***
## fe.123 0.52520655 0.02053143 25.58062 < 2e-16 ***
## fe.124 0.64677946 0.01899700 34.04640 < 2e-16 ***
## fe.125 0.70745885 0.01899072 37.25287 < 2e-16 ***
## fe.126 0.74050153 0.02037904 36.33643 < 2e-16 ***
## fe.127 0.75699649 0.01969678 38.43249 < 2e-16 ***
## fe.128 0.74561411 0.02028314 36.76029 < 2e-16 ***
## fe.129 0.76356595 0.01987532 38.41779 < 2e-16 ***
## fe.130 0.74716447 0.01972862 37.87211 < 2e-16 ***
## fe.131 0.70897756 0.01920377 36.91867 < 2e-16 ***
## fe.132 0.61333371 0.01847301 33.20161 < 2e-16 ***
## fe.133 0.49802008 0.01901317 26.19343 < 2e-16 ***
## fe.134 0.52475971 0.01789866 29.31838 < 2e-16 ***
## fe.135 0.57090833 0.01838795 31.04795 < 2e-16 ***
## fe.136 0.58133085 0.01845842 31.49407 < 2e-16 ***
## fe.137 0.61702259 0.01831299 33.69317 < 2e-16 ***
## fe.138 0.61982246 0.01842202 33.64574 < 2e-16 ***
## fe.139 0.54661787 0.01870485 29.22333 < 2e-16 ***
## spei_lag 0.00362411 0.00442644 0.81874 0.41294
## ndvi_lag -0.00247898 0.00318826 -0.77753 0.43685
## smoist_lag 0.95780169 0.00415816 230.34290 < 2e-16 ***
## sp_spei 0.00281464 0.00582238 0.48342 0.62880
## sp_ndvi -0.08618798 0.00442928 -19.45871 < 2e-16 ***
## sp_smoist -0.13823932 0.00612426 -22.57241 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.303305 on 48644 degrees of freedom
## Number of observations: 48789 Degrees of Freedom: 48644
## SSR: 4474.94813 MSE: 0.091994 Root MSE: 0.303305
## Multiple R-Squared: 0.906683 Adjusted R-Squared: 0.906407Predictive performance
To evaluate the predictive capabilities of our model, we fit it on the first 300 months of data. We also fit 2 simpler models (VAR and SAR) with which we will compare ours. More details on this are in Section 4.2.
# Prepare the dataset containing only the first 300 months
df_short <-
prepare_spvar_df(test_df, test_df_spei, test_df_smoist, test_df_ndvi, T =
301)
spvar_fit_short <-
systemfit(
c(
eqspei = spei ~ 0 + . -spei -ndvi -smoist,
eqndvi = ndvi ~ 0 + . -spei -ndvi -smoist,
eqsmoist = smoist ~ 0 + . -spei -ndvi -smoist
),
data = df_short,
method = "SUR"
)
var_fit_short <-
systemfit(
c(
eqspei = spei ~ 0 + . -spei -ndvi -smoist -sp_spei -sp_smoist -sp_ndvi,
eqndvi = ndvi ~ 0 + . -spei -ndvi -smoist -sp_spei -sp_smoist -sp_ndvi,
eqsmoist = smoist ~ 0 + . -spei -ndvi -smoist -sp_spei -sp_smoist -sp_ndvi
),
data = df_short,
method = "SUR"
)
sar_fit_short <-
systemfit(
c(
eqspei = spei ~ 0 + . -spei -ndvi -smoist -sp_smoist -sp_ndvi -ndvi_lag -smoist_lag,
eqndvi = ndvi ~ 0 + . -spei -ndvi -smoist -sp_smoist -sp_spei -spei_lag -smoist_lag,
eqsmoist = smoist ~ 0 + . -spei -ndvi -smoist -sp_spei -sp_ndvi -ndvi_lag -spei_lag
),
data = df_short,
method = "SUR"
)We now look at the adjusted R^2 of the three models for the SPEI equation.
print(paste("The adjusted R Squared of the SpVAR is:",adj_r2(spvar_fit_short$eq[[1]])))## [1] "The adjusted R Squared of the SpVAR is: 0.82613029802549"print(paste("The adjusted R Squared of the VAR is:",adj_r2(var_fit_short$eq[[1]])))## [1] "The adjusted R Squared of the VAR is: 0.822843111752807"print(paste("The adjusted R Squared of the SAR is:",adj_r2(sar_fit_short$eq[[1]])))## [1] "The adjusted R Squared of the SAR is: 0.823165740366314"We subset the data as needed for making predictions and set all the missing values to 0 so that they don’t impact our predictions.
data_pred_short <- spvar_df[((299 * 139) + 1):dim(spvar_df)[1],]
data_pred_short[is.na(data_pred_short)] = 0We now look at predictive performances of the three models using the \(SMAPE\) and the \(RMSE\). We start by comparing the error rates using the three models for predictions up to a year in advance
for (i in c(1, 3, 6, 12)) {
print(paste(i, "steps ahead:"))
spvar_err <-
compute_pred_steps(102 - i, i, spvar_fit_short, data_pred_short)
print(paste("The SMAPE and the RMSE are as follows for the SpVAR:", spvar_err[[2]],spvar_err[[1]]))
var_err <-
compute_pred_steps(102 - i, i, var_fit_short, data_pred_short, i+1)
print(paste("The SMAPE and the RMSE are as follows for the VAR:", var_err[[2]],var_err[[1]] ))
sar_err <-
compute_pred_steps(102 - i, i, sar_fit_short, data_pred_short, i)
print(paste("The SMAPE and the RMSE are as follows for the SAR:", sar_err[[2]], sar_err[[1]]))
}## [1] "1 steps ahead:"
## [1] "The SMAPE and the RMSE are as follows for the SpVAR: 0.519449953439214 0.354568699970673"
## [1] "The SMAPE and the RMSE are as follows for the VAR: 0.544368037695199 0.373058356467052"
## [1] "The SMAPE and the RMSE are as follows for the SAR: 0.53543895506312 0.365162180946665"
## [1] "3 steps ahead:"
## [1] "The SMAPE and the RMSE are as follows for the SpVAR: 0.841867634373804 0.632879764211257"
## [1] "The SMAPE and the RMSE are as follows for the VAR: 0.907297560555501 0.716630402536862"
## [1] "The SMAPE and the RMSE are as follows for the SAR: 0.900163896789106 0.711343337616749"
## [1] "6 steps ahead:"
## [1] "The SMAPE and the RMSE are as follows for the SpVAR: 1.13616085581854 0.865707314827308"
## [1] "The SMAPE and the RMSE are as follows for the VAR: 1.21334638604964 1.08265638615537"
## [1] "The SMAPE and the RMSE are as follows for the SAR: 1.20819845037974 1.09194854260833"
## [1] "12 steps ahead:"
## [1] "The SMAPE and the RMSE are as follows for the SpVAR: 1.43219968494042 1.04396822410734"
## [1] "The SMAPE and the RMSE are as follows for the VAR: 1.46570276462238 1.67487168383349"
## [1] "The SMAPE and the RMSE are as follows for the SAR: 1.47879993729122 1.87327111660701"We check if the differences are significant using the Diebold-Mariano test.
for (i in c(1, 3, 6, 12)) {
spvar_res <-
compute_pred_steps(102 - i,
i,
spvar_fit_short,
data_pred_short,
value = "res")
sar_res <-
compute_pred_steps(102 - i,
i,
sar_fit_short,
data_pred_short,
i,
value = "res")
var_res <-
compute_pred_steps(102 - i,
i,
var_fit_short,
data_pred_short,
i + 1,
value = "res")
spvar_res <- spvar_res[2:((102 - i) * 44)]
sar_res <- sar_res[2:((102 - i) * 44)]
var_res <- var_res[2:((102 - i) * 44)]
print(dm.test(spvar_res, sar_res, h = i))
print(dm.test(spvar_res, var_res, h = i))
}##
## Diebold-Mariano Test
##
## data: spvar_ressar_res
## DM = -2.4919, Forecast horizon = 1, Loss function power = 2, p-value =
## 0.01274
## alternative hypothesis: two.sided
##
##
## Diebold-Mariano Test
##
## data: spvar_resvar_res
## DM = -5.0626, Forecast horizon = 1, Loss function power = 2, p-value =
## 4.304e-07
## alternative hypothesis: two.sided
##
##
## Diebold-Mariano Test
##
## data: spvar_ressar_res
## DM = -9.128, Forecast horizon = 3, Loss function power = 2, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##
##
## Diebold-Mariano Test
##
## data: spvar_resvar_res
## DM = -9.8962, Forecast horizon = 3, Loss function power = 2, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##
##
## Diebold-Mariano Test
##
## data: spvar_ressar_res
## DM = -14.681, Forecast horizon = 6, Loss function power = 2, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##
##
## Diebold-Mariano Test
##
## data: spvar_resvar_res
## DM = -13.679, Forecast horizon = 6, Loss function power = 2, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##
##
## Diebold-Mariano Test
##
## data: spvar_ressar_res
## DM = -19.269, Forecast horizon = 12, Loss function power = 2, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##
##
## Diebold-Mariano Test
##
## data: spvar_resvar_res
## DM = -19.414, Forecast horizon = 12, Loss function power = 2, p-value <
## 2.2e-16
## alternative hypothesis: two.sidedWe now fit the models on a shorter data-set (containing only the first 150 months) to see the performance on longer term predictions (60 and 120-months ahead).
# Prepare the dataset containing only the first 150 months
df_super_short <-
prepare_spvar_df(test_df, test_df_spei, test_df_smoist, test_df_ndvi, T =
151)
spvar_fit_super_short <-
systemfit(
c(
eqspei = spei ~ 0 + . -spei -ndvi -smoist,
eqndvi = ndvi ~ 0 + . -spei -ndvi -smoist,
eqsmoist = smoist ~ 0 + . -spei -ndvi -smoist
),
data = df_super_short,
method = "SUR"
)
var_fit_super_short <-
systemfit(
c(
eqspei = spei ~ 0 + . -spei -ndvi -smoist -sp_spei -sp_ndvi -sp_smoist,
eqndvi = ndvi ~ 0 + . -spei -ndvi -smoist -sp_spei -sp_ndvi -sp_smoist,
eqsmoist = smoist ~ 0 + . -spei -ndvi -smoist -sp_spei -sp_ndvi -sp_smoist
),
data = df_super_short,
method = "SUR"
)
sar_fit_super_short <-
systemfit(
c(
eqspei = spei ~ 0 + . -spei -ndvi -smoist -sp_smoist -sp_ndvi -ndvi_lag -smoist_lag,
eqndvi = ndvi ~ 0 + . -spei -ndvi -smoist -sp_smoist -sp_spei -spei_lag -smoist_lag,
eqsmoist = smoist ~ 0 + . -spei -ndvi -smoist -sp_spei -sp_ndvi -ndvi_lag -spei_lag
),
data = df_super_short,
method = "SUR"
)We again subset the data as needed and set all the missing values to 0 so that they don’t impact our predictions.
data_pred_super_short <- spvar_df[((149 * 139) + 1):dim(spvar_df)[1],]
data_pred_super_short[is.na(data_pred_super_short)] = 0Finally, we check the performances on these longer-term predictions and check if the differences between models are significant with the Diebold-Mariano test.
for (i in c(60, 120)) {
print(paste(i, "steps ahead:"))
spvar_err <-
compute_pred_steps(252 - i, i, spvar_fit_super_short, data_pred_super_short)
print(paste(
"The SMAPE and RMSE are as follows for the SpVar",
spvar_err[[2]],
spvar_err[[1]]
))
var_err <-
compute_pred_steps(252 - i, i, var_fit_super_short, data_pred_super_short, i+1)
print(paste(
"The SMAPE and RMSE are as follows for the VAR",
var_err[[2]],
var_err[[1]]
))
sar_err <-
compute_pred_steps(252 - i, i, sar_fit_super_short, data_pred_super_short, i)
print(paste(
"The SMAPE and RMSE are as follows for the SAR",
sar_err[[2]],
sar_err[[1]]
))
spvar_res <-
compute_pred_steps(252 - i,
i,
spvar_fit_super_short,
data_pred_super_short,
value = "res")
sar_res <-
compute_pred_steps(252 - i,
i,
sar_fit_super_short,
data_pred_super_short,
i,
value = "res")
var_res <-
compute_pred_steps(252 - i,
i,
var_fit_super_short,
data_pred_super_short,
i + 1,
value = "res")
spvar_res <- spvar_res[2:((252 - i) * 139)]
sar_res <- sar_res[2:((252 - i) * 139)]
var_res <- var_res[2:((252 - i) * 139)]
print(dm.test(spvar_res, sar_res), h = i)
print(dm.test(spvar_res, var_res), h = i)
}## [1] "60 steps ahead:"
## [1] "The SMAPE and RMSE are as follows for the SpVar 1.49133291703587 1.03689778569715"
## [1] "The SMAPE and RMSE are as follows for the VAR 1.52336136792023 1.9454493372135"
## [1] "The SMAPE and RMSE are as follows for the SAR 1.55726484587559 2.67271385831612"
##
## Diebold-Mariano Test
##
## data: spvar_ressar_res
## DM = -50.84, Forecast horizon = 1, Loss function power = 2, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##
##
## Diebold-Mariano Test
##
## data: spvar_resvar_res
## DM = -49.32, Forecast horizon = 1, Loss function power = 2, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##
## [1] "120 steps ahead:"
## [1] "The SMAPE and RMSE are as follows for the SpVar 1.5289622005954 1.08035932682959"
## [1] "The SMAPE and RMSE are as follows for the VAR 1.53116762201427 2.33478155223426"
## [1] "The SMAPE and RMSE are as follows for the SAR 1.57837456210087 3.32292953764326"
##
## Diebold-Mariano Test
##
## data: spvar_ressar_res
## DM = -41.619, Forecast horizon = 1, Loss function power = 2, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##
##
## Diebold-Mariano Test
##
## data: spvar_resvar_res
## DM = -41.788, Forecast horizon = 1, Loss function power = 2, p-value <
## 2.2e-16
## alternative hypothesis: two.sidedImpulse Response Functions
We now look at impulse response functions starting from same region IRFs.
reg_2 <- 106
t_irf <- 100
quantiles <- bootstrap(t_irf, reg_2, spvar_fit, spvar_df, 1000)
multiplot_irf_sameregion(spvar_fit, spvar_df, reg_2, 100, quantiles)
We now look at the IRF in neighbouring regions when shocking a region having transitional climate.
plot <- multiplot_irf_all(spvar_fit, 100, reg_2, plot_all = FALSE)
map <- multiplot_irf_map(spvar_fit, index_pairs, test_df, 50, reg_2, 20)
map[[1]] <-
map[[1]] + theme(plot.margin = unit(c(
b = 0.3,
l = 0,
t = 0.3,
r = 0
), "cm"))
map[[2]] <-
map[[2]] + theme(plot.margin = unit(c(
b = 0.3,
l = 0,
t = 0.3,
r = 0
), "cm"))
map[[3]] <-
map[[3]] + theme(plot.margin = unit(c(
b = 0.3,
l = 0,
t = 0.3,
r = 0
), "cm"))
grid.arrange(
plot[[1]],
map[[1]],
plot[[2]],
map[[2]],
plot[[3]],
map[[3]],
nrow = 3,
ncol = 2,
widths = c(1, 1.1)
)
Finally, we plot the IRFs for neighbouring regions when shocking regions with alpine and Mediterranean climate.
reg_1 <- 19
reg_3 <- 137
plot_1 <- multiplot_irf_all(spvar_fit, 75, reg_1, plot_all = FALSE)
plot_2 <- multiplot_irf_all(spvar_fit, 75, reg_3, plot_all = FALSE)
grid.arrange(plot_1[[1]], plot_2[[1]], plot_1[[2]], plot_2[[2]],
plot_1[[3]], plot_2[[3]],
ncol = 2)
Other plots
The following plots the mean monthly precipitations and the average PET for each region.
# Prepare a dataframe for plotting
plot_df <-
create_plot_df(1,
reg_latitudes ,
reg_long_min,
reg_long_max,
lats,
lons,
SPEI,
NDVI,
SMoist)
# Compute the mean precipitation and PET and store them in the dataframe
for (i in 1:139) {
lon_ind <- as.numeric(strsplit(index_pairs[i], ",")[[1]][2])
lat_ind <- as.numeric(strsplit(index_pairs[i], ",")[[1]][1])
plot_df[plot_df$lat_ind == lat_ind &
plot_df$lon_ind == lon_ind, 5] <-
mean(PET[lon_ind, lat_ind, 6:402])
plot_df[plot_df$lat_ind == lat_ind &
plot_df$lon_ind == lon_ind, 6] <-
mean(PRE[lon_ind, lat_ind, 6:402])
}
#Duplicate the dataframe
plot_df_2 <- plot_df
#Rename a column for plotting
colnames(plot_df_2)[5] <- "weights"
colnames(plot_df)[6] <- "weights"
plot_1 <-
create_plot(plot_df_2,
var = "Mean PET", scale_fixed = TRUE)
plot_2 <-
create_plot(plot_df,
var = "Mean PRE", scale_fixed = TRUE)
grid.arrange(plot_1, plot_2, ncol = 2)
The following plots the evolution of PET and PRE (monthly precipitations) in certain regions (one for each climatic area) to highlight the differences between them.
# Prepare a dataframe with the data
test_df_pet <- create_test_df(reg_latitudes,
reg_long_min,
reg_long_max,
lats,
lons,
PET,
PRE,
SPEI)
# Rename relevant columns
names(test_df_pet)[6] <- "PET"
names(test_df_pet)[7] <- "PRE"
plot_1_pre <-
ggplot(data = test_df_pet[test_df_pet$index == index_pairs[reg_1] |
test_df_pet$index == index_pairs[reg_2] |
test_df_pet$index == index_pairs[reg_3], ]) +
geom_smooth(aes(x = month, y = PET, color = index)) +
guides(colour = "none") +
scale_color_manual(
name = "Climate",
labels = c("Mediterranean", "Transitional", "Alpine"),
values = c("274,384", "271,384", "262,392")
) + labs(x = "Months")
plot_2_pre <-
ggplot(data = test_df_pet[test_df_pet$index == index_pairs[reg_1] |
test_df_pet$index == index_pairs[reg_2] |
test_df_pet$index == index_pairs[reg_3],]) +
geom_smooth(aes(x = month, y = PRE, color = index)) + scale_color_manual(
name = "Climate",
labels = c("Mediterranean", "Transitional", "Alpine"),
values = c("274,384", "271,384", "262,392")
) + labs(x = "Months")
grid.arrange(plot_1_pre, plot_2_pre, ncol = 2, widths = c(0.9, 1.3))
The following plots the region under consideration (Continental Italy) and how we divide it to exemplify some concepts of spatial analysis together with the weighting scheme we use.
# Prepare a plot showcasing the region under consideration
plot_df<-create_plot_df(1,reg_latitudes, reg_long_min, reg_long_max, lats, lons, SPEI, SMoist, NDVI)
reg_plot <- qmplot(
y = lat_plot,
x = lon_plot,
data = plot_df,
geom = "blank",
maptype = "toner-background",
zoom = 7
) +
geom_rect(
aes(
ymin = lat_min,
ymax = lat_max,
xmin = lon_min,
xmax = lon_max,
alpha = 0.01,
color = "white",
fill = "grey"
)
) + scale_alpha_continuous(guide = "none") +
scale_color_discrete(guide = "none", type = "white") +
scale_fill_discrete(guide = "none", type = "grey") +
annotate(
"text",
x = 12.75,
y = 46.75,
label = "B",
size = 4,
color = "black"
) +
annotate(
"text",
x = 11.75,
y = 45.25,
label = "A",
size = 4,
color = "black"
) +
annotate(
"text",
x = 15.75,
y = 40.75,
label = "C",
size = 4,
color = "black"
)
# Prepare a plot showcasing the weighting scheme
plot_weights <- create_plot_df_weights(weights[reg_2,])
weight_plot <- create_plot(plot_weights, var = "Weights")
reg_plot <-
reg_plot + theme(plot.margin = unit(c(0, 0.75, 0, 0), "cm"))
grid.arrange(reg_plot,
weight_plot,
ncol = 2,
widths = c(0.87, 1))
The following plots the distribution of \(\lambda\) under the first and tenth unit root test
lambda_first_test <- tibble(x = lambda_mc[1:10000, ])
lambda_tenth_test <- tibble(x = lambda_mc[90001:100000, ])
hist1 <-
ggplot(data = lambda_first_test) + geom_histogram(
aes(x = x),
binwidth = 0.0005,
color = "black",
fill = "black"
) + labs(x = "Lambda in test 1", y = "Count") +
geom_vline(aes(xintercept = quantile(x, probs = 0.05)), color = "red") +
scale_y_continuous(limits = c(0, 570))
hist2 <-
ggplot(data = lambda_tenth_test) + geom_histogram(
aes(x = x),
binwidth = 0.0005,
color = "black",
fill = "black"
) + labs(x = "Lambda in test 10", y = NULL) +
scale_y_continuous(limits = c(0, 570), guide = "none") +
geom_vline(aes(xintercept = quantile(x, probs = 0.05)), color = "red")
grid.arrange(hist1, hist2, ncol = 2, widths = c(1.2, 1))
The following plots the number of droughts of each type in the 402 months under consideration.
plots_droughts <-
long_term_prediction_plot(
test_df_spei[1:402,], plot="initial"
)
The following plots the predicted number of droughts in the second 201 periods of the period under consideration (below) against the actual number of them in the data (above). Note that we use a model fit only on the first 201 months.
# Prepare the data
test_df_201<- test_df[((139 * 199) + 1):(139 * 201), ]
test_df_201$month <- c(rep(1, times = 139), rep(2, times = 139))
test_df_201_spei <- test_df_spei[200:201, ]
test_df_201_ndvi <- test_df_ndvi[200:201, ]
test_df_201_smoist <- test_df_smoist[200:201, ]
# Fit a model on the first 201 periods
spvar_df_201 <- spvar_df[1:(201 * 139), ]
spvar_fit_201 <-
systemfit(
c(
eqspei = spei ~ 0 + . - spei - ndvi - smoist,
eqndvi = ndvi ~ 0 + . - spei - ndvi - smoist,
eqsmoist = smoist ~ 0 + . - spei - ndvi - smoist
),
data = spvar_df_201,
method = "SUR"
)
# Create the plot
plots_droughts_201 <-
long_term_prediction_plot(
test_df_spei[201:402, ],
test_df_201,
test_df_201_spei,
test_df_201_smoist,
test_df_201_ndvi,
t_pred = 201,
spvar_fit = spvar_fit_201
)
The following plots the predicted number of droughts in the next 402 periods (the 402 months after those covered by the dataset, i.e. starting from January 2014).
# Create a table with the data of the last 2 years
test_df_last_year <- create_test_df(
reg_latitudes,
reg_long_min,
reg_long_max,
lats,
lons,
SPEI_last_year,
SMoist_last_year,
NDVI_last_year
)
# Subset the data
test_df_last_year <- test_df_last_year[1:278, ]
#Create a separate df for each variable
test_df_last_year_spei <-
data.frame(split(test_df_last_year$SPEI, as.factor(test_df_last_year$index)))
test_df_last_year_ndvi <-
data.frame(split(test_df_last_year$NDVI, as.factor(test_df_last_year$index)))
test_df_last_year_smoist <-
data.frame(split(
test_df_last_year$SMoist,
as.factor(test_df_last_year$index)
))
plots_droughts_future <-
long_term_prediction_plot(
test_df_spei[201:402, ],
test_df_last_year,
test_df_last_year_spei,
test_df_last_year_smoist,
test_df_last_year_ndvi,
t_pred = 402,
spvar_fit = spvar_fit,
plot = "future"
)